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'?.  ABSTRACT  (Maximum  200  words)  1 

This  report  documents  the  development  and  testing  of  nowcasting  methods  for  satellite 
imagery.  Techniques  for  identifying,  tracking,  and  forecasting  features  of  interest  through  contours  of  I 
brightness  temperature  in  geostationary  IR  imagery  were  evaluated  in  a  study  covering  over  200  j 

features  from  all  four  seasons.  Results  indicate  that  because  contours  frequently  split,  merge,  or  j 

dissipate,  contour  based  techniques  are  successful  in  only  a  minority  of  cases.  The  most  skillful  ■ 

contour  extrapolation  technique  is  movable  persistence,  which  simply  extrapolates  the  position  of  the  j 

contour  based  on  the  past  movement  of  its  centroid  position.  Its  verification  scores  are  significantly  ; 

better  than  persistence  for  forecasts  from  1/2  to  3  hours  in  advance. 

As  an  alternative  to  the  contour-based  schemes,  a  cross-correlation  technique  for  predicting 
short-term  changes  in  pixel  grayshade  values  (brightness  temperatures)  was  developed  and  tested. 
Pixel-based  verification  of  the  cross-correlation  forecasts  show  an  improvement  over  persistence  out  to 
2  1/2  hours  (the  longest  time  period  studied  for  this  technique).  Contour  forecasts  derived  from  the  ' 

cross-correlation  are  significantly  more  skillful  than  those  of  the  best  contour  extrapolation  technique  ! 

(movable  persistence)  out  to  1  hour,  and  no  worse  out  to  2  1/2  hours.  j 
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One  of  the  weather  forecast  tools  wanted  by  the  U.S.  Air  Force  is  one  to 
help  make  more  effective  nowcasts  of  cloud  and  precipitation  features.  In  the 
future,  the  weather  forecaster  will  have  more  data  available  than  can  be 
synthesized  in  a  short  amount  of  time,  so  schemes  which  simplify  the  job  of 
forecasting  are  desirable.  Another  desirable  characteristic  for  such  a  scheme 
is  the  ability  to  work  with  very  little  data,  since  many  Air  Force  operations 
are  in  data-denied  areas  or  areas  where  active  sensors  cannot  be  operated. 
The  Air  Force  Phillips  Laboratory  has  been  actively  engaged  in  developing 
and  testing  prototype  nowcasting  schemes  meeting  these  criteria.  Once 
developed  and  adequately  tested,  a  prototype  scheme  could  eventually  be 
transferred  to  an  Air  Force  meteorological  workstation  such  as  the 
Automated  Weather  Distribution  System,  or  AWDS. 

The  focus  of  previous  work  at  the  Phillips  Laboratory  for  nowcasting 
clouds  and  precipitation  has  been  on  contour  extrapolation  techniques  for 
radar  and  satellite  data.  The  basic  idea  of  the  contour  extrapolation  schemes 
is  as  follows:  the  cloud  elements  with  significant  weather  can  be  identified  as 
areas  bounded  by  some  significant  contour;  for  example,  the  precipitating 
area  of  a  thunderstorm  may  be  defined  by  a  cloud  top  temperature  of  colder 
than  -50  C.  With  a  time  series  of  bounding  contours  around  the  -50  C  area, 
an  intelligent  extrapolation  of  the  future  size,  shape,  and  movement  may  be 
estimated  (Figure  1).  A  number  of  forecast  schemes  have  been  developed  or 
tailored  at  the  Phillips  Laboratory.  They  are:  (1)  the  whole  contour  technique 
(Bohne  et  al.  1988)  (Bianco  and  Huang  1990);  (2)  the  segmentation  method 
(ibid);  (3)  the  adaptive  exponential  smoothing  technique,  also  called  the 


Kawas  technique  (Kawa?  and  Chen  1989),  and  (4)  movable  persistence. 
Section  2  of  this  report  will  elaborate  more  on  these  schemes. 


Prior  testing  of  these  schemes  had  been  limited  by  the  cumbersome 
contour  selection  process,  which  had  not  been  automated  and  involved 
contouring  a  satellite  image,  producing  paper  plots  of  the  contours,  searching 
the  paper  plots  for  trackable  features,  etc.  Some  problems  were  quickly 
discovered:  first,  this  method  was  very  slow.  Second,  many  of  the  desired 
features  did  not  remain  coherent  over  a  large  number  of  frames,  i.e.,  one 
feature  split  into  many,  two  features  merged  into  one,  a  feature  totally 
dissipated.  None  of  the  contouring  schemes  were  designed  to  handle  such 
cases.  Because  of  these  limitations,  a  comparison  of  forecast  skill  among  the 
different  schemes  proved  inconclusive  due  to  small  sample  size  (Heideman  et 
al.  1990;  Ruggiero  et  al.  1991). 


CONTOURS  CONTOUR 


Figure  1.  Figure  illustrating  how  initial  frames  (time  t-At  and  t,  on  left) 
from  a  time  series  of  analyzed  contours  can  be  used  to  generate 
forecast  contours  (for  time  t+At,  on  right). 


2 


As  a  result,  the  Phillips  Lab  initiated  a  contract  to  adequately  test  the 
contouring  schemes  and  to  recommend  one  winner  for  technology  transition 
or  further  development.  This  report  documents  the  work  performed  under 
this  contract.  As  agreed  with  the  Phillips  Laborato**y,  we  proposed  in  the 
contract  to: 

(1)  Develop  a  prototype  display  and  archive  software  to  allow  a  user  to 
loop  up  consecutive  frames  of  remapped  satellite  imagery  or  display  each 
separately.  Using  a  mouse,  the  user  could  then  select  a  pixel  having  the 
brightness  temperature  desired  to  use  for  a  contour,  and  display/overlay 
contours  for  each  of  the  frames  of  satellite  imagery.  Last,  this  software  would 
archive  the  contours  for  use  by  the  forecasting  and  verification  programs. 

(2)  Examine  possible  enhancements  to  the  existing  algorithms  which 
would  broaden  their  generality.  Ways  to  deal  with  the  splitting  and  merging 
of  contours,  to  forecast  multiple  cloud  layers,  and  to  deal  with  boundary 
intersections  would  be  studied. 

(3)  Study  replacement  algorithms  which  may  be  more  generally 
applicable.  These  algorithms  would  be  abandoned  or  recommended  as 
candidates  for  further  development  and  testing,  depending  on  the  results  of 
the  test  with  the  contouring  algorithms. 

(4)  Test  each  algorithm  over  a  broad  range  of  cases  (50  for  each 
season)  and  verify  algorithm  skill  using  the  standard  skill  scores  such  as 
probability  of  detection,  false  alarm  rate,  critical  success  index,  and  root 
mean  square  error.  Selection  of  the  cases  to  be  tested  was  to  be  done  on  the 
basis  of  meteorological  interest,  not  on  what  the  algorithm  could  handle. 
Rephrased,  to  simulate  the  process  in  a  real  forecasting  situation,  the  areas 
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to  be  selected  had  dynamic  or  interesting  weather.  This  was  done  to  simulate 
the  situations  in  which  a  weather  forecaster  was  likely  to  seek  forecast  advice 
from  an  automated  algorithm.  Thus,  if  any  or  all  of  the  schemes  did  not 
work,  this  would  have  been  noted.  The  results  would  indicate  the  percent  of 
the  time  each  scheme  was  able  to  produce  a  usable  forecast.  Further,  if 
generalities  could  be  made  about  regimes  where  an  algorithm  is  strong  or 
weak,  these  would  be  made. 

One  modification  that  was  made  shortly  after  we  began  the  contract 
was  to  abandon  the  search  for  algorithmic  improvements  which  could  handle 
splitting  and  merging  and  instead  focus  on  the  development  of  a  new  forecast 
scheme  based  on  cross-correlations.  This  was  done  for  two  reasons.  First, 
after  some  preliminary  research,  we  saw  no  easy  way  to  handle  the  splitting 
and  merging  of  features,  and  attempting  work  on  this  problem  could  have 
slowed  down  the  evaluation  process.  Further,  with  four  candidate  forecast 
techniques,  four  general  solutions  might  have  been  required,  but  the  short 
timeline  of  the  project  may  have  dictated  our  improving  only  one  of  the  four 
techniques,  thus  possibly  biasing  the  evaluation  to  the  one  selected  for 
improvement.  The  second  major  reason  was  that  the  technique  of  cross¬ 
correlations  (Bellon  et  al.  1992)  avoided  many  of  the  pitfalls  of  the  contouring 
techniques,  producing  full  forecast  IR  satellite  images  rather  than  contour 
forecasts.  Thus,  time  allotted  for  tasks  (2)  and  (3)  above  were  mostly  spent  in 
the  development  of  a  new,  prototype  cross-correlations  scheme. 

This  report,  then,  documents  the  resulting  contract  work  performed 
over  the  last  year  on  satellite  nowcasting.  Some  of  the  work,  namely,  the 
testing  of  the  cross-correlations  technique  on  synoptic- scale  features,  has 
already  been  described  (Hamill  and  Nehrkom  1993).  In  this  report,  we  will 
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first  give  a  quick  overview  of  the  forecast  schemes  (Section  2),  and  an  in- 
depth  summary  of  the  evaluation  process  and  the  results  (Section  3).  The 
last  part,  Section  4,  provides  conclusions. 

2.  Description  of  Contour  Selection  and  Contour  Forecast 
Algorithms 

This  section  describes  the  new  contour  selection  and  cross-correlations 
forecast  software  developed  this  year  and  reviews  the  existing  contour 
forecast  schemes. 

2.1  Contour  Selection  Software 

A  way  of  selecting  trackable  contour  features  was  obviously  needed 
before  forecasts  could  be  run  or  verified.  As  mentioned  before,  the  pre¬ 
existing  process  for  selecting  cases  for  testing  the  algorithms  was 
cumbersome,  involving  a  lengthy  remapping  process,  production  of  paper 
plots  of  contours,  and  other  time-consuming  steps.  With  new  software  tools 
available  (specifically,  an  image  processing  language  called  PVWAVE)  it  was 
decided  the  first  part  of  the  time  should  be  allotted  to  developing  new  front- 
end  software.  This  software  would  be  an  integrated  package  to  interactively 
remap  satellite  data,  allow  the  user  to  examine  the  remapped  data  and  to 
select  features  for  examination,  contour  them,  and  store  the  contour  and  pixel 
information  if  desired  for  later  input  to  the  forecast  algorithms.  The 
components  of  this  system  are  summarized  below. 

2.1.1  Remapping  of  Geostationary  Data 

The  forecast  domain  selected  from  which  cases  were  chosen  is  shown  in 
Figure  2.  It  is  a  Lambert-conformal  map  projection  true  at  30  and  60 
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degrees  N,  with  the  domain  centered  at  41.36  degrees  N  latitude,  and  88.68 
degrees  W  longitude  (Marseilles,  IL).  The  new  software  remaps  to  a  256*256 
grid,  with  an  effective  resolution  of  approximately  11  km  (later  the  grid  is 
expanded  to  a  512*512  grid  through  pixel  replication).  Contour  selection  and 
evaluation  were  judged  to  be  easier  and  more  accurate  in  a  Lambert 
projection  than  in  the  projection  of  the  raw  GOES  data,  where  the  earth's 
curvature  was  apparent. 

Previously,  remapping  GOES  data  to  the  Lambert  projection  was  done 
in  batch  mode  because  of  the  extended  length  of  the  job.  The  remapping 
process  now  has  been  speeded  up  so  it  can  be  done  interactively.  The  old 
remapping  software  was  designed  to  perform  a  coordinate  conversion  for 
every  pixel  on  the  Lambert  grid,  a  computationally  intensive  process.  The 
algorithm  would  start  with  pixel  (1,1),  determine  the  location  in  the  raw 
GOES  image,  and  remap  the  pixel  value,  and  proceed  to  (2,1),  (3,1),  and  so  on. 
Because  the  coordinate  conversion  can  be  considered  linear  over  small 
sections  of  the  forecast  domain,  we  decided  that  only  a  select  few  points 
needed  to  be  fully  coordinate  converted,  and  the  rest  of  the  conversion  could 
be  done  by  interpolation.  For  example,  say  that  we  know  pixel  (1,1) 
corresponds  to  pixel  (a,b)  in  the  forecast  image,  and  (3,1)  corresponds  to  (c,d). 
Then,  assuming  linearity,  we  can  calculate  the  coordinate  of  (2,1)  as  ( (a+c)/2  , 
(b+d)/2  ).  In  our  algorithm,  we  actually  perform  the  coordinate  conversion 
every  15  pixels  in  each  direction,  and  determine  the  rest  through  bilinear 
interpolation.  The  resulting  remapped  image  is  indistinguishable  from  the 
same  image  remapped  with  a  full  coordinate  conversion,  and  the  time  savings 
is  tremendous;  on  the  VAX  3600  workstation,  remapping  six  frames  of 
imagery  takes  2-3  minutes,  whereas  the  old  method  took  30  minutes. 
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2.1.2  Contour  Level  Selection. 

After  the  image  has  been  remapped,  the  analyst  has  a  set  of  tools 
available  for  inspecting  the  image  and  selecting  an  interesting  subimage  and 
a  level  for  contouring.  Some  of  the  features  are: 

(1)  Looping  capabilities. 

(2)  Ability  to  examine  a  single  frame  of  the  image. 

(3)  Enhancement  tables. 

(a)  GOES  ZA  and  ZB  tables  (Clark  1983), 

(b)  A  set  of  15  false  color  enhancement  tables, 

(c)  Interactive  generation  of  new  enhancement  tables. 

(4)  Storage  of  remapped  imagery  for  quicker  reprocessing  later. 

The  user  also  has  a  choice  of  how  a  contour  level  can  be  selected.  The 
user  can:  (1)  Enter  the  pixel  brightness  value  to  use  as  a  contour,  or  (2)  Click 
on  a  pixel  on  the  screen  and  use  the  value  of  the  pixel  at  that  location.  The 
user  also  can  select  the  size  of  the  subdomain  around  the  feature,  inside 
which  the  contouring  will  take  place.  Figure  3  illustrates  one  sample 
remapped  frame  of  satellite  imagery  that  might  be  examined  for  selection  of  a 
trackable  contour  feature. 

2.1.3  Contour  Generation  /  Selection 

The  new  package  is  also  set  up  to  contour  the  image  and  to  display  the 
resulting  contours  overlaid  on  the  original  satellite  image  loop  (Figure  4). 
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Since  there  may  be  more  than  one  contour  in  a  given  subframe,  if  the  user 
decides  that  one  of  the  features  is  trackable  (i.e.,  contiguous  over  many 
frames),  then  the  software  will  prompt  the  user  to  click  a  mouse  button  inside 
the  feature  of  interest.  The  software  will  automatically  delete  any  other 
contours,  and  once  again  display  a  satellite  loop  with  the  one  selected  contour 
feature  overlaid  (Figure  5). 

2.1.4  Storage  of  Contour  and  Pixel  Information. 

The  last  step  performed  is  the  storage  of  either  contour  information 
(for  the  contour  forecast  techniques)  or  actual  satellite  subframes  (for  the 
cross-correlation  scheme).  The  contour  information  is  stored  in  a  form  that 
is  usable  as  input  by  all  the  contouring  schemes,  and  the  satellite  subframes 
can  be  read  into  the  correlation  scheme  directly. 

2.2  Forecast  Algorithm  Description 

The  evaluation  of  the  cloud  forecasts  was  to  be  done  for  six  techniques; 
the  first  five  were  all  contouring  techniques,  and  the  last,  the  cross¬ 
correlations,  produces  pixel  forecast  images,  from  which  contours  can  be 
extracted  for  comparison  with  the  other  schemes.  A  brief  explanation  of  each 
of  the  schemes  is  given  below. 

2.2.1  Whole  Contour  Technique 

The  whole  contour  technique  is  more  completely  described  in  Bianco 
and  Huang  (1990).  The  premise  of  this  scheme  is  that  a  contour  can  by 
sampled  at  regular  points  along  the  contour,  and  described  by  a  periodic  set 
of  x-  and  y-coordinates.  A  graph  of  the  x-  and  y-components  shows  they  are 
periodic;  thus,  they  can  be  decomposed  by  Fourier  analysis  into  a  spectral 
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representation  of  phases  and  amplitudes  (Figure  6).  The  number  of  waves 
describing  the  contour  can  be  specified  by  the  user;  the  more  waves,  the  more 
accurate  the  initial  description  of  the  shape  of  the  contour. 

In  addition  to  phases  and  amplitudes,  the  centroid  positions  and  aspect 
ratios  (length  /  height)  are  also  calculated.  All  the  calculated  parameters, 
including  the  time  tendencies  for  each  Fourier  component,  are  then  linearly 
extrapolated  in  time,  and  a  forecast  contour  is  reconstructed  through  an 
inverse  transform.  The  forecast  shape  and  size  are  scaled  to  fit  the 
extrapolated  aspect  ratio  and  area.  The  position  of  the  contour  is  determined 
from  the  extrapolated  centroid  position,  with  the  centroid  position  defined  as 
the  average  of  the  contour  boundary  values. 

The  testing  of  this  scheme  allowed  the  user  to  vary  the  number  of 
previous  observations  used  to  initialize  the  model.  For  two  previous 
observations,  only  a  simple  linear  extrapolation  was  possible.  With  three  or 
more  observations,  this  linear  extrapolation  estimate  was  improved  through  a 
least-squares  regression.  The  user  was  also  able  to  vary  the  number  of 
Fourier  components  used  in  the  analysis  (from  1  to  5)  and  the  forecast 
method,  which  allowed  linear  extrapolation  of  the  Fourier  amplitude  and 
phases,  or  more  complex  extrapolations  including  relationships  between 
adjacent  Fourier  components. 

2.2.2  Segmentation  Technique 

The  segmentation  technique  is  also  described  in  Bianco  and  Huang 
(1990),  where  it  is  referred  to  as  the  "angle  displacement  method."  The 
abbreviated  description  given  in  Ruggiero  et  al.  (1991)  is  used  here. 
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The  segmentation  technique  defines  a  contour  by  its  centroid  position 
and  the  length  of  lines  radiating  out  from  the  centroid  to  the  contour 
boundary  (Figure  7).  The  lines  are  at  equally  spaced  angles  around  the 
centroid.  The  angle  separating  the  line  segments  is  called  the  displacement 
(alpha)  angle  and  is  determined  by  the  user.  The  segmentation  technique 
uses  a  linear  least  squares  fit  to  forecast  future  line  lengths  and  centroid 
position.  The  contour  is  reconstructed  by  connecting  the  endpoints  of  the 
lines  radiating  out  from  the  forecast  centroid. 

Displacement  angles  of  5,  10  ,15,  and  20  degrees  were  tested. 

2.2.3  Adaptive  Exponential  Smoothing  (Kawas)  Technique 

The  Kawas  technique  is  very  similar  to  the  segmentation  technique. 

It  is  described  in  Kawas  and  Chen  (1989).  This  technique  uses  the  centroid 
position  and  line  lengths  radiating  from  the  centroid  to  describe  the  contour. 
For  Kawas,  however,  the  displacement  angle  is  fixed  at  22.5  degrees. 
Extrapolation  is  done  by  an  adaptive  exponential  smoothing  procedure  which 
is  a  weighted  least  squares  fitting  procedure  with  an  exponentially 
discounted  weighting  factor. 

The  user  was  required  to  set  a  smoothing  coefficient  to  determine  how 
sensitive  the  method  is  to  rapid  changes  in  the  observations.  A  small  value  of 
the  smoothing  coefficient  puts  more  emphasis  on  the  older  observations, 
while  a  large  value  makes  the  extrapolation  more  reliant  on  newer 
observations. 
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2.2.4  Movable  Persistence 

This  is  the  simplest  of  the  techniques.  Shape  and  size  of  the  feature 
are  persisted,  but  the  centroid  position  is  extrapolated  from  the  latest 
position  using  a  linear  least-squares  regression  centroid  displacement  in 
previous  frames. 

After  much  of  the  forecast  work  was  done,  it  was  determined  that  the 
methodology  for  defining  centroids  could  be  enhanced.  Rather  than  defining 
a  centroid  as  the  average  of  the  contour  locations,  the  centroid  was  defined  as 
the  center  of  mass.  In  this  way,  long,  thin  features  which  were  not 
contributing  a  large  amount  of  volume  to  the  overall  contour  were  not  heavily 
weighted  in  the  centroid  position. 

As  will  be  explained  later,  verification  scores  for  movable  persistence 
were  found  to  be  erroneous  after  completion  of  the  forecasts.  Corrected 
movable  persistence  scores  were  recomputed  for  both  methods  of  defining  the 
centroid  position.  For  logistical  reasons,  scores  were  only  computed  for 
forecasts  initialized  from  the  first  two  frames  of  the  12-frame  sequence  of 
images  used  in  the  case  selection. 

2.2.5  Persistence 

Straight  persistence  of  the  last  analyzed  contour  is  also  included  for 
comparison  in  the  verification  statistics.  As  is  discussed  in  the  results 
section,  verification  scores  for  persistence  were  found  to  be  erroneous  for  the 
majority  of  cases  because  of  a  coding  error.  Corrected  persistence  scores 
were  recomputed  along  with  those  for  movable  persistence,  again  only  for 
forecasts  initialized  from  the  first  two  frames  of  the  12-frame  sequence  of 
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images  used  in  the  case  selection.  Comparison  with  the  other  techniques  is 
thus  restricted  to  this  subsample. 

2.2.6  Cross  -  Correlations  Technique 

This  new  technique  was  introduced  to  work  around  the  limitations  of 
the  contour-based  forecast  techniques.  As  discussed  in  the  introduction,  all 
straight  contour  forecast  techniques  have  the  disadvantage  of  being  unable  to 
handle  cases  where  forecast  features  split,  merge,  or  dissipate.  Additionally, 
an  implicit  assumption  of  these  schemes  is  that  the  important  information 
from  a  satellite  scene  can  be  synthesized  into  a  series  of  contours.  This 
results  in  a  binary  cloud  forecast;  a  point  is  either  inside  the  contour  or 
outside;  actual  brightness  temperature  or  rainfall  rate  is  not  forecast. 

Bellon  et  al.  (1992)  describe  a  possible  approach  to  cloud  forecasting, 
one  which  can  implicitly  handle  splitting  and  merging  and  which  does  result 
in  a  full  forecast  image  of  IR  pixel  brightness  values.  Judging  from  the 
apparent  success  of  this  approach,  we  decided  to  attempt  developing  such  a 
scheme  at  the  Phillips  Laboratory.  A  description  of  the  technique  and  its 
application  to  the  synoptic  scale  was  described  in  Atch.  1  (Hamill  and 
Nehrkom  1993).  We  review  some  of  the  background  from  that  paper  below, 
and  point  out  how  the  technique  is  applied  differently  to  the  small-scale 
features  used  in  the  contour  extrapolation  studies. 

At  the  heart  of  this  prototype  cloud  forecast  technique  is  the  derivation 
of  displacement  vectors  (i.e.,  "winds")  through  a  cross  correlations  analysis. 
First,  a  subset  of  pixels  from  the  first  image  in  a  satellite  loop  is  chosen.  For 
purpose  of  illustration  in  Figure  8,  this  subset  is  8*8;  in  our  actual  prototype 
scheme,  the  subset  chosen  is  15*15.  Next,  based  on  the  maximum  possible 
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wind  displacement  in  half  an  hour,  a  search  radius  is  chosen,  and  identically 
sized  subsets  of  pixels  from  the  second  image  with  centers  inside  the  search 
radius  are  correlated  against  the  subset  from  the  first  image.  If  the  scene  has 
a  uniform  cloud  layer,  then  the  advective  velocity  for  this  layer  will  be  defined 
by  the  vector  from  the  center  of  the  subset  in  Image  1  to  the  center  of  the 
subset  with  the  maximum  correlation  in  Image  2.  This  is  illustrated  in 
Figure  8,  and  a  sample  plot  of  correlation  coefficients  and  a  derived 
displacement  vector  is  shown  in  Figure  9.  In  our  technique,  a  cross 
correlations  analysis  is  repeated  at  a  gridded  subset  of  points  throughout  the 
domain,  and  wind  displacement  vectors  for  all  points  are  then  derived 
through  an  objective  analysis.  A  backwards-in-time  trajectory  forecast  is 
then  used  to  produce  the  cloud  forecasts. 

There  are  a  number  of  differences  in  how  the  technique  has  been 
tailored  to  handle  the  small-scale  features  rather  than  the  synoptic-scale. 
First,  the  GOES  data  used  in  the  synoptic-scale  forecast  experiment  was  on  a 
256*256  grid,  with  an  approximate  resolution  of  11  km.  For  the  selection  of 
smaller  features  to  track  and  forecast,  this  domain  was  enlarged  to  512*512, 
halving  the  effective  resolution  of  the  satellite  data.  Second,  the  quality 
control  step  was  modified  slightly  to  throw  out  any  correlation  vector  with  a 
maximum  correlation  below  0.75.  Subjectively,  it  was  found  that  this  was  an 
effective  compromise  between  eliminating  bad  displacement  vectors  and 
having  few  wind  vectors.  Another  modification  was  the  use  of  three 
iterations  in  the  objective  analysis  instead  of  two,  with  the  radius  of  influence 
interactively  determined  from  the  size  of  the  feature  being  forecast;  the 
smaller  the  feature,  the  smaller  the  radii  of  influence,  and  the  greater  the 
detail.  Otherwise,  the  forecast  algorithm  is  identical  to  the  synoptic-scale 
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algorithm  described  in  the  attached  conference  paper.  Figure  10  shows  the 
raw,  quality-controlled,  and  objectively  analyzed  displacement  vectors  used  to 
derive  the  subsequent  trajectory  forecast  for  the  satellite  picture  shown  in 
Figure  5. 

In  addition  to  the  pixel-by-pixel  verification  of  forecast  satellite  images 
discussed  in  the  appendix,  we  also  derived  computed  contours  from  the  cross¬ 
correlation  forecasts.  These  contour  forecasts  were  then  compared  to  the 
contour  extrapolation  schemes. 

3.  Evaluation  of  forecast  techniques 

3.1.  Approach 

Testing  of  the  contour  extrapolation  techniques  posed  some  special 
problems,  since  the  testing  environment  is  quite  different  from  the  eventual 
operational  environment  in  which  the  candidate  techniques  might  ultimately 
be  used.  We  attempted  to  make  the  process  of  selecting  the  contour  levels 
and  resulting  features  as  objective  as  possible.  Procedures  used  were 
designed  to  guard  against  the  tendency  to  select  only  those  features  that  were 
trackable  and  resulted  in  usable  cases  for  the  verification  statistics.  A  set  of 
12  consecutive  half-hourly  satellite  images  were  loaded  for  each  case  selection 
process.  Only  the  first  6  images  were  used  to  select  one  or  more  features  for 
tracking  and  forecasting.  The  last  6  images  were  considered  to  be  future 
verification  images  not  available  to  the  user.  The  forecasts  required  at  least 
3  successive  non-splitting,  non-merging,  and  non-dissipating  observations  -  at 
least  2  observations  from  the  first  6  images  for  the  extrapolation  and  at  least 
1  observation  from  the  last  6  images  for  verification.  If  successive 
observations  were  available  only  in  the  first  6  images,  the  selection  was 
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considered  a  failure  because  of  a  lack  of  at  least  1  verification  observation. 
Though  the  last  6  images  were  treated  as  verification  images  in  the  case 
selection  process,  observations  available  from  these  images  were  also  used  to 
produce  forecasts  which  used  more  than  6  observations  in  the  extrapolation. 

The  percentage  of  case  selection  failures  for  each  season  were:  fall  - 
18%;  winter  -  26%;  spring  -  25%;  summer  -  24%.  The  fall  failure  rate  is  lower 
because  the  above  case  selection  criteria  were  not  utilized  for  this  season 
until  after  65%  of  the  cases  had  been  selected.  The  percentage  of  splitting, 
merging,  and  dissipating  contours  which  caused  the  above  failures  were:  fall  - 
42%,  33%,  25%;  winter  -  47%,  26%,  27%;  spring  -  35%,  47%,  18%;  and 
summer  -  62%,  25%,  13%.  It  should  be  noted  that  even  if  a  case  selection  was 
successful  by  our  definition,  the  selected  feature  was  trackable  through  all 
available  subsequent  image  frames  in  only  a  minority  of  cases.  The 
percentage  of  splitting,  merging,  and  dissipating  contours  which  prevented 
the  use  of  observations  out  to  the  12th  image  for  the  successful  case  selections 
were:  fall  (92%)  -  36%,  28%,  28%;  winter  (67%>  27%,  20%,  20%;  spring  (69%)- 
21%,  23%,  25%;  and  summer  (75%)  -  27%,  15%,  33%.  On  average,  the 
selected  features  were  identifiable  in  slightly  over  half  of  all  available  frames. 
The  number  of  frames  for  each  season  were:  fall  -  6.4;  winter  -  7.6;  spring  - 
7.7;  summer  -  7.2. 

There  are  a  large  number  of  possible  combinations  of  parameters  for  which 
verification  scores  can  be  computed.  We  distinguish  between  two  basic  types 
of  parameters:  characteristics  of  the  cases  (size  and  speed  of  the  features,  and 
the  number  of  previous  observations  available  for  forecasting),  and 
extrapolation  technique  parameters.  We  chose  a  subset  of  all  possible 
combinations,  based  on  the  twin  objectives  of  this  part  of  the  study:  (1) 
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determination  of  optimal  parameter  settings  for  contour  extrapolation 
techniques,  and  (2)  intercomparison  of  the  four  extrapolation  techniques,  and 
comparison  with  other  techniques  (cross-correlation  and  persistence).  For 
this  control  s  t  of  cases,  sensitivity  tests  were  performed  by  systematically 
varying  the  adjustable  parameters  of  each  forecast  scheme.  We  repeated 
these  sensitivity  tests  for  other  subsamples  of  cases  to  determine  whether  the 
results  depend  on  the  attributes  of  the  cases  (e.g.,  whether  different 
parameter  values  give  the  best  results  for  slow  and  fast  moving  features). 
Based  on  the  results  of  these  tests,  optimal  parameter  value(s)  for  each 
scheme  were  determined.  Once  the  optimal  parameter  values  had  been 
determined,  they  were  held  fixed  for  the  remainder  of  all  forecast  method 
intercomparisons  (but  we  allowed  different  optimal  values  for  different 
seasons,  if  appropriate).  Forecast  method  intercomparisons  concentrated  on 
the  control  set  of  cases,  but  other  samples  (i.e.,  slow  and  fast  features,  small 
and  large  features)  were  also  used  to  determine  whether  the  ranking  of  the 
methods  depends  on  the  case  characteristics. 

The  main  verification  score  used  for  comparison  was  the  Critical 
Success  Index  (CSI),  although  four  others  (Probability  of  Detection  -POD, 
False  Alarm  Rate  -  FAR,  Bias,  and  Root  mean  square  error  of  centroid 
position  -  RMSE)  were  als  *  considered.  As  illustrated  in  Figure  11,  the  scores 
can  be  defined  as  follows: 

CSI  =  A  /  (A+B+C) 

FAR  =  C  /  (C+A) 


POD  =  A  /  (B+A) 


Bias  = 


A  +  C 
A  +  B 


where  A  represents  the  cloud  pixels  correctly  forecast,  C  represents  false 
alarms  (forecast,  no  verification),  and  B  the  misses  (verification,  no  forecast). 
Scores  for  CSI,  FAR,  and  POD  range  from  0.0  to  1.0;  a  high  CSI  and  POD  is 
desirable,  as  is  a  low  FAR. 

In  the  following,  we  discuss  the  stratification  of  the  statistics  by  case 
parameters,  and  describe  the  details  of  the  forecast  method  sensitivity  and 
intercomparison  tests. 

3.2.  Stratification  of  cases 


Cases  are  characterized  by  four  parameters:  observed  speed  and  size  of 
the  feature,  the  number  of  previous  frames  available  for  forecasting,  and  the 
forecast  increment  or  lead  time.  In  this  definition  of  "case”,  each  forecast  is 
considered  a  separate  case,  whereas  in  the  database  each  case  number  refers 
to  all  forecasts  for  a  particular  feature.  We  formed  three  categories  each  for 
the  size  and  speed  of  the  features  (small;  medium;  large  and  slow;  moderate; 
fast,  respectively).  Selection  criteria  for  these  categories  were  chosen  such 
that  between  V3  and  V2  of  all  cases  ended  up  in  tne  middle  categories.  The 
number  of  previous  frames  available  for  forecasting  ranges  from  2  to  11,  and 
the  forecast  increment  ranges  from  1  to  6  frames  (V2  to  3  hours).  The  control 
set  of  cases  (CONTROL)  consisted  of  all  cases  with: 

size:  medium  or  large 
speed:  slow,  moderate,  or  fast 

number  of  previous  frames:  any  number  between  2  and  11. 
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That  is,  all  cases  except  those  with  small  features  are  included  in 
CONTROL.  Small  features  were  excluded  to  eliminate  features  that  are  less 
likely  to  be  of  meteorological  significance.  For  an  examination  of  how  much 
the  results  depend  on  the  case  characteristics,  we  also  considered  samples  for 
which  one  of  three  case  parameters  is  varied  from  CONTROL.  The  total 
number  of  samples  is  summarized  in  the  following  table: 


Table  1:  Case  parameter  values  of  samples  (sets  of  cases)  used  for  forecast 
method  testing  and  verification. 


Name 

Size 

Speed 

Previous  Frames 

CONTROL 

medium,large 

(CONTROL) 

slow,  moderate, 

fast 

(CONTROL) 

2-11 

(CONTROL) 

SMALL 

small 

CONTROL 

CONTROL 

LARGE 

large 

CONTROL 

CONTROL 

SLOW 

CONTROL 

slow 

CONTROL 

FAST 

CONTROL 

fast 

CONTROL 

PRE2 

CONTROL 

CONTROL 

PRE4 

CONTROL 

CONTROL 

4 

PRE6 

CONTROL 

CONTROL 

6 

18 


The  above  table  represents  a  total  of  8  samples.  For  each  sample,  all 
available  lead  times  were  considered.  Because  of  the  limited  lifetime  of 
features  and  limitations  on  the  number  of  frames  available  for  case  selection, 
not  all  combinations  of  number  of  previous  frames  and  forecast  increments 
are  present  for  each  feature.  Thus,  sample  size  and  composition  are  not  the 
same  for  different  forecast  increments  -  this  is  acceptable  since  the  crucial 
intercomparison  is  between  different  methods  and  forecast  parameter 
settings,  at  each  forecast  increment,  not  between  different  forecast 
increments. 


The  distribution  of  size  and  speed  of  the  features  was  examined  in 
detail  for  the  fall  season.  Sizes  for  the  53  features  span  more  than  2  orders  of 
magnitude  (52  to  9100  km2) .  When  plotted  on  a  semi-logarithmic  plot,  the 
distribution  is  bimodal,  with  one  peak  at  fairly  small  sizes  (ca.  200  km2),  and 
one  just  above  the  median  (ca.  500  km2).  Some  percentiles  are  shown  in 
Table  2: 

Table  2:  Distribution  of  feature  sizes  for  the  fall  season 


rank 

size 

(km2) 

rank 

size 

(km2) 

V6: 

9 

95 

5/q: 

45 

1990 

V4: 

13 

152 

2/4: 

41 

1224 

V3: 

18 

207 

V3: 

35 

661 

V2: 

27 

433 

(median) 
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We  used  cutoff  values  corresponding  to  the  approximate  V4  -  3 U 
thresholds,  and  roughly  centered  around  the  median  (on  a  logarithmic  scale): 


small  (11  features):  features  no  larger  than  150  km2 
large  (13  features):  features  larger  than  1200  km2 

Use  of  these  thresholds  removes  only  the  very  smallest  features  from  the 
CONTROL  sample,  leaving  a  sufficient  number  to  allow  the  further 
stratifications  indicated  in  Table  1. 

Speeds  have  a  unimodal  distribution  which  is  skewed  towards  slow 
speeds,  and  a  long  tail  at  high  speeds.  Some  percentiles  of  the  53  fall 
features  are  shown  in  Table  3: 


Table  3:  Distribution  of  translation  speeds  for  the  fall  season 


rank 

speed 

(km/h) 

rank 

speed 

(km/h) 

V6: 

9 

55 

5/e' 

45 

135 

V4: 

13 

60 

3/4: 

41 

113 

V3: 

18 

68 

2/3: 

35 

105 

V2: 

27 

88 

(median) 

When  all  the  small  features  are  excluded,  the  distribution  with  speed 
for  the  remaining  cases  is  modified  as  follows: 

Total  number  of  cases:  42 
Median:  Rank=21.5  Speed=89  km/h 

To  assure  sample  sizes  of  at  least  10  for  the  SLOW  and  FAST  samples, 
we  used  approximate  V3  -  2/3  cutoffs  centered  around  the  median  value: 

slow  (14  features):  features  with  speeds  no  faster  than  70  km/h 
last  (12  features):  features  with  speeds  faster  than  110  km/h 

For  the  winter  season,  the  size  and  speed  distributions  of  the  features 
were  similar.  The  total  number  of  52  cases  were  subdivided  into  small,  large 
and  medium  features  according  to  the  following  cutoff  values: 

small  (11  features):  features  no  larger  than  250  km2 
large  (11  features):  features  larger  than  1000  km2 

After  removal  of  the  small  features,  the  remaining  41  features  were 
further  subdivided  into  slow,  moderate,  and  fast  speeds  according  to  the 
following  cutoff  values: 

slow  (15  features):  features  with  speeds  no  faster  than  80  km/h 
fast  (12  features):  features  with  speeds  faster  than  120  km/h 

For  spring  and  summer,  scores  were  only  computed  for  the  control 
sample  (examination  of  the  results  by  subsample  revealed  no  significant 
differences  between  the  fall  and  winter  seasons).  We  used  a  cutoff  of  250  km2 
(225  km2)  for  the  small  cases  for  spring  (summer),  eliminating  10  (13) 
features  from  the  total  sample  of  52  (51),  leaving  42  (38)  features  in  the 
control  sample. 
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3.3.  Sensitivity  tests 


3.3.1.  Whole  Contour  Method: 

There  are  three  separate  tunable  parameters  for  this  technique:  the 
number  of  Fourier  modes  used  for  prediction  (frequency),  the  type  of  forecast 
method,  and  the  order  of  extrapolation.  However,  previous  experimentation 
has  already  established  linear  extrapolation  (order=l)  as  the  optimal  choice. 
The  remaining  two  parameters  have  the  following  possible  values: 

frequency:  1-5 
method:  2,  3,  or  4 

The  verification  scores  for  the  CONTROL  sample  for  fall  are  displayed 
in  Figure  12  for  all  values  of  frequency  for  method  2.  The  corresponding 
scores  for  method  3  and  4  are  shown  in  Figures  13  and  14.  The  sample  sizes 
for  the  CONTROL  sample  range  from  499  (increment  1),  337  (increment  2), 
207  (increment  3),  114  (increment  4),  60  (increment  5),  to  22  (increment  6).  A 
striking  feature  of  the  plots  of  the  Critical  Success  Index  (CSI)  is  the 
insensitivity  of  the  results  to  frequency.  Close  examination  reveals,  however, 
that  for  all  methods,  for  all  increments,  best  CSIs  are  obtained  for 
frequency=l  (with  one  exception,  where  frequency=2).  For  those  frequency 
values,  CSI  is  completely  insensitive  to  method.  This  conclusion  also  applies 
to  the  probability  of  detection  (POD),  false  alarm  rate  (FAR),  and  the  bias. 
The  rms  error  of  the  centroid  position  (RMSE)  is  insensitive  to  frequency  and 
method. 

It  is  somewhat  surprising  that  using  only  one  frequency  mode,  i.e., 
approximating  the  shape  of  any  cloud  feature  selected  for  tracking  by  an 
ellipse,  results  in  better  forecasts  than  retaining  and  predicting  higher 
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wavenumbers.  It  indicates  that  there  is  no  skill,  at  least  using  the 
approaches  tried  with  the  whole  contour  method,  in  predicting  any  but  the 
grossest  scales  of  cloud  features. 

Figures  15  and  16  show  the  CSI  and  RMSE  for  the  LARGE  and 
SMALL  subsamples,  respectively  (because  there  is  essentially  no  sensitivity 
to  the  forecast  method,  only  results  for  method  no.  2  are  shown  in  these  and 
all  following  figures).  Results  agree  qualitatively  with  those  for  the 
CONTROL  sample,  as  do  those  for  the  SLOW  and  FAST  samples  (Figures  17 
and  18)  and  the  samples  with  2 , 4,  or  6  previous  frames  (PRE2,  PRE4,  and 
PRE6),  corresponding  to  Figures  19,  20,  and  21.  Sample  sizes  for  the 
subsamples  are  substantially  smaller  than  for  CONTROL,  ranging  from  100  - 
200  at  increment  1  to  5  - 15  at  increment  6;  sample  sizes  for  the  PRE4  and 
PRE6  subsamples  are  even  smaller,  with  no  cases  for  the  longer  increments. 

For  the  purposes  of  forecast  method  intercomparison,  we  chose  an 
optimum  parameter  setting  based  on  the  results  described  above  of 
frequency:*!.,  and  method  =  2  (although  any  one  of  the  three  methods  (2,  3,  or 
4)  would  give  the  same  results). 

Comparing  skill  between  subsamples  gives  some  interesting  results. 
CSI  values  increase  with  the  size  of  the  feature  (they  are  higher  for  LARGE 
than  for  CONTROL,  and  smaller  for  SMALL  than  for  CONTROL).  The 
RMSE  are  much  smaller  for  small  features,  but  in  light  of  the  extremely 
small  CSI  (0.2)  the  overall  forecast  skill  for  small  features  must  still  be 
judged  poor;  LARGE  and  CONTROL  RMSEs  are  roughly  comparable,  with  a 
tendency  towards  smaller  errors  for  large  features.  Both  measures  of  skill 
improve  as  the  speed  of  the  feature  decreases  (values  for  the  SLOW  sample 
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are  better  than  those  for  CONTROL,  which  in  turn  are  better  than  those  for 
the  FAST  sample).  There  is  no  clear  dependence  of  the  CSI  or  RMSE  skills 
on  the  number  of  previous  frames.  This  is  not  surprising,  since  only  the 
previous  2  frames  are  used  for  forecasting  in  our  linear  extrapolation 
technique. 

Results  for  the  other  three  seasons  are  shown  in  Figures  22,  23,  and  24 
for  the  control  sample.  The  lack  of  sensitivity  to  the  number  of  modes  is 
evident  in  these  seasons  as  well.  CSI  scores  tend  to  be  slightly  higher  in  fall 
and  winter  than  spring  and  summer,  particularly  at  longer  forecast 
increments. 

3.3.2.  Segmentation  Method 

The  only  tunable  parameter  for  this  technique  is  the  segmentation 
angle,  which  was  varied  from  5°  to  20°  in  increments  of  5°  (a  total  of  4 
values). 

Skill  scores  for  the  CONTROL  sample  for  fall  are  shown  in  Figure  25. 
There  is  very  little  sensitivity  with  respect  to  the  segmentation  angle  in 
either  the  CSI  or  the  RMSE,  with  some  suggestion  that  neither  the  minimum 
nor  the  maximum  value  (5°,  20°)  is  optimal.  The  POD  generally  decreases 
with  increasing  segmentation  angle,  whereas  the  FAR  and  the  bias  is  largest 
at  the  smallest  angle,  but  nearly  constant  otherwise  for  most  forecast 
increments.  Thus,  the  higher  values  of  POD  at  smaller  segmentation  angles 
come  at  the  expense  of  overforecasting.  The  RMSE  is  quite  insensitive  to  the 
segmentation  angle,  aside  from  increment  6,  where  minimum  values  are  for 
segmentation  angles  of  10°  and  15°. 
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The  lack  of  sensitivity  to  the  segmentation  angle  indicates  that  for  this 
technique,  as  was  found  for  the  whole  contour  technique,  there  is  no  skill  in 
predicting  small  scale  features  of  the  selected  contours.  Visual  inspection  of 
segmentation  forecasts  revealed  that  smaller  segmentation  angles  not  only 
fail  to  improve  the  verification  scores,  but  also  result  in  highly  unrealistic 
contour  shapes.  Since  adjacent  segment  lengths  are  forecast  independently 
from  one  another,  the  forecasted  contours  can  contain  long,  thin  protrusions, 
and  bear  no  resemblance  to  meteorological  features.  Since  the  area  enclosed 
by  these  protrusions  is  small,  however,  the  skill  scores  are  not  degraded 
significantly  for  these  unrealistic  contour  forecasts.  This  results  in 
deceptively  good  scores  in  some  cases. 

The  results  for  the  SLOW  and  FAST  subsamples  (Figures  26  and  27) 
are  qualitatively  similar  to  those  for  the  CONTROL  sample.  The  LARGE 
sample  results  (Figure  28)  also  confirm  the  results  of  the  CONTROL  sample, 
whereas  the  SMALL  sample  (Figure  29)  indicates  better  results  for  larger 
angles  (15  or  20  degrees). 

We  chose  a  value  of  15°  as  the  optimal  parameter  setting  for  the 
purpose  of  forecast  method  intercomparison  based  on  these  results. 

The  differences  in  skill  scores  between  the  subsamples  mirror  those 
discussed  previously  for  the  whole  contour  method.  In  particular,  skill 
increases  with  increasing  size  and  decreasing  speed  of  the  feature,  and  shows 
no  systematic  dependence  on  the  number  of  previous  frames  available  for 
forecasting  (Figures  30  -  32;  recall  that  the  segmentation  method  also  only 
uses  the  previous  2  frames  for  a  linear  extrapolation). 
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Results  for  the  other  three  seasons  are  shown  in  Figures  33,  34,  and  35 
for  the  control  sample.  The  lack  of  sensitivity  to  the  segmentation  angle  is 
evident  in  these  seasons  as  well.  CSI  scores  tend  to  be  slightly  higher  in  fall 
and  winter  than  spring  and  summer,  particularly  at  longer  forecast 
increments,  as  was  seen  in  the  whole  contour  technique. 

3.3.3.  Kawas  Method: 

The  only  tunable  parameter  for  this.  'chmque  is  the  smoothing 
coefficient,  which  was  varied  from  0.10  to  O.bO  in  increments  of  .05. 

For  the  CONTROL  sample  for  fall  (Figure  36),  the  CSI,  POD,  and  FAR 
all  indicate  increasing  skill  with  increasing  smoothing  coefficient  for  forecast 
increments  less  than  3.  At  longer  forecast  increments,  optimal  values  are 
less  than  at  shorter  increments  ( .5  for  increments  1  and  2,  .35  and  .4  for 
increments  3-6).  The  bias  is  smallest  for  smoothing  coefficients  of  .4  to  .5  for 
short  increments  (1-3),  whereas  it  increases  monotonically  with  smoothing 
coefficient  for  longer  increments.  The  RMSE  indicates  smallest  errors  for 
smoothing  coefficients  in  the  range  .25  to  .35,  for  all  forecast  increments. 

The  sensitivity  of  the  CSI  and  RMSE  with  respect  to  the  smoothing 
coefficient  are  similar  for  the  SLOW  and  FAST  subsample  (Figures  37  and 
38).  However,  optimum  values  are  slightly  larger  for  fast  features,  at  least  in 
terms  of  the  RMSE  (.3  to  .4,  vs.  .25  to  .35).  The  SLOW  subsample  shows  just 
the  opposite,  with  optimum  values  slightly  lower  (.3  to  .45  for  CSI,  .15  to  .25 
for  RMSE  for  increments  >  1). 

The  SMALL  and  LARGE  subsample  (Figures  39  and  40)  indicate 
similar  opposing  sensitivities:  CSI  for  small  features  increase  with  smoothing 
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coefficient  (optimum  values  >  .5),  for  large  features  optimum  values  are  .3  to 
.45.  The  RMSE  for  large  features  indicates  optimum  values  of  .2  to  .25, 
whereas  in  the  SMALL  subsample  optimum  values  range  from  .15  for 
increment  1  to  .45  for  increment  6. 

There  are  no  striking  differences  in  sensitivity  for  the  pre2,  pre4,  or 
pre6  subsamples  (Figures  41  -  43). 

In  contrast  to  the  other  two  techniques,  skill  scores  show  a  much 
greater  sensitivity  to  the  tunable  parameter  of  the  technique.  In  theory,  this 
might  present  an  opportunity  for  optimization  of  the  technique;  in  practice, 
however,  it  will  complicate  its  implementation.  Choosing  an  overall  optimum 
value  for  the  forecast  method  intercomparison  is  difficult  for  the  Kawas 
technique,  since  different  optimum  values  exists  for  the  different  increments, 
subsamples,  and  measures  of  skill.  However,  it  is  possible  to  choose  one 
value  that  results  in  scores  that  are  not  too  different  from  the  optimum  score 
for  any  increment,  for  either  CSI  or  RMSE.  We  selected  an  optimum  value  of 
.40  based  on  the  results  of  the  CONTROL  sample.  This  value  is  also  close  to 
the  optimum  values  for  the  LARGE  and  SLOW  subsamples.  For  the  SMALL 
and  FAST  subsamples  it  is  decidedly  suboptimal,  particularly  at  short 
forecast  increments;  however,  since  the  CSI  skill  for  these  subsamples  is  low 
for  any  value  of  the  smoothing  coefficient,  this  is  not  a  serious  problem. 

As  was  the  case  in  the  other  two  techniques,  CSI  scores  are  generally 
higher  for  larger  features,  and  RMSE  errors  tend  to  increase  with  feature 
size.  For  smoothing  coefficients  at  the  respective  optimum  values,  CSI  and 
RMSE  scores  tend  to  improve  as  the  speed  decreases. 
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Comparison  between  the  PRE2,  PRE4,  and  PRE6  subsamples  indicates 
no  clear  advantage  of  having  a  longer  time  history  on  which  to  base  the 
prediction. 

Results  for  the  other  three  seasons  are  shown  in  Figures  44,  45,  and  46 
for  the  control  sample.  The  sensitivity  to  the  smoothing  coefficient  is  similar 
in  these  seasons.  Because  in  the  winter  sample  slightly  better  scores  are 
obtained  for  larger  values  of  the  smoothing  coefficient,  we  used  an  optimum 
value  of  .45  for  winter  for  subsequent  method  intercomparisons.  In  spring 
and  summer,  we  used  the  same  value  as  in  fall  (.40). 

3.4.  Forecast  method  intercomparison 

3.4.1  Intercomparison  of  three  contour  extrapolation  techniques 

Figure  47  shows  a  direct  comparison  of  the  skill  scores  for  the 
CONTROL  sample  for  fall  of  the  contour  extrapolation  methods,  using  the 
optimum  values  selected  above  (whole  contour:  method=2,  frequency=l; 
segmentation:  angle=15°;  Kawas:  smoothing  coefficient  =  0.4).  Also  shown 
are  the  sample  sizes  for  the  different  forecast  increments  and  techniques 
(sample  sizes  for  the  different  techniques  are  not  exactly  the  same  because 
the  extrapolation  algorithms  failed  for  some  cases  for  one  of  the  techniques, 
but  not  the  others;  these  differences  in  sample  size  are  small  enough, 
however,  to  avoid  influencing  the  results).  It  can  be  seen  that  the  segmenta¬ 
tion  technique  is  superior  overall  in  terms  of  CSI  and  POD;  its  FAR  is  slightly 
inferior  to  the  whole  contour  method  at  larger  increments,  and  its  bias  is 
worse  than  that  of  both  other  methods.  The  Kawas  is  inferior  to  both 
methods  in  terms  of  CSI,  POD,  and  FAR;  its  bias,  however,  is  the  smallest  of 
any  of  the  forecast  methods.  Thus,  the  segmentation  method  tends  to 
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overforecast  more  than  whole  contour,  but  the  improved  POD  more  than 
makes  up  for  it  in  the  CSI,  and  the  Kawas  method  is  overall  inferior  by  any 
measure  except  bias.  In  terms  of  the  RMSE,  the  segmentation  method  is 
better  at  most  increments  than  the  whole  contour  RMSE,  and  the  Kawas 
method  is  by  far  worse  than  the  other  two  methods. 

The  relative  ranking  in  terms  of  CSI  and  RMSE  of  the  three  methods 
is  the  same  in  the  CONTROL,  SLOW,  FAST,  LARGE  and  SMALL 
subsamples  (Figures  48,  49,  50  and  51):  the  segmentation  is  best,  followed  by 
the  whole  contour  method,  with  Kawas  a  distant  third  (particularly  in  terms 
of  RMSE). 

Results  for  the  other  three  seasons  confirm  the  conclusion  that  the 
segmentation  method  is  the  best  of  the  three  contour  extrapolation  methods 
considered  so  far  (Figures  52,  53,  54).  The  relative  ranking  of  the  three 
techniques  for  all  the  measures  of  skill  is  the  same  for  fall,  winter,  and 
spring.  In  summer,  the  Kawas  method  is  second  to  segmentation,  ahead  of 
the  whole  contour  method,  mainly  because  its  FAR  in  that  season  is  better 
when  compared  to  the  other  techniques. 

3.4.2.  Significance  Testing  and  comparison  with  other  techniques 

We  performed  a  statistical  analysis  on  the  forecast  method 
intercomparison  to  assess  the  statistical  significance  of  the  results.  We 
separately  considered  each  score  for  each  forecast  increment.  An  analysis  of 
variance  (ANOVA;  see  Neter  and  Wassermann,  1974,  hereafter  referred  to  as 
NW)  was  performed  on  the  raw  verification  scores,  using  the  Splus  statistical 
package  (Chambers  and  Hastie,  1992).  Because  there  is  substantial 
variability  of  the  forecast  score  from  one  case  to  the  next  for  each  of  the 
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methods,  we  considered  the  forecast  method  a  treatment  factor,  and  the  case 
number  a  blocking  factor.  Different  forecasts  produced  for  the  same  case 
number  were  grouped  together  under  each  case  number.  The  treatment 
factor  (forecast  method)  is  a  fixed  effect,  since  the  methods  were  chosen 
beforehand  for  their  intrinsic  interest,  whereas  the  blocking  factor  must  be 
considered  a  random  effect,  since  the  cases  are  a  sample  from  a  larger 
population.  For  this  analysis  we  eliminated  all  cases  that  did  not  have  scores 
for  all  forecast  methods  available,  so  that  a  complete  block  design  resulted 
(all  methods  were  tested  for  each  of  the  blocks,  i.e.  cases).  With  this 
experimental  design,  tests  for  the  presence  of  block  and  treatment  effects  are 
performed  using  a  two-way  ANOVA  (NW,  p.  474  and  p.  730).  In  the  presence 
of  treatment  effects,  the  confidence  interval  for  all  pairwise  comparisons 
between  treatment  means  (i.e.,  mean  scores  for  the  different  forecast 
methods)  can  be  computed  using  the  Tu1  ,  confidence  intervals  (NW,  p.  474 
and  p.  730).  Results  from  this  analysis  are  presented  below.  Assumptions 
underlying  these  significance  tests  are  that  the  residuals  are  normally 
distributed  random  variables  (residuals  are  the  deviations  of  the  actual 
scores  from  the  means  for  the  method  and  case  number),  and  that  the  mean 
scores  for  the  cases  are  normally  distributed  random  variables. 

We  first  applied  the  ANOVA  to  a  comparison  of  the  three  contour 
extrapolation  techniques  with  persistence  and  movable  persistence.  As  was 
discussed  in  sections  2.2.4  and  2.2.5,  scores  for  persistence  and  movable 
persistence  were  only  available  for  forecasts  initialized  from  the  first  two 
frames  of  the  12-frame  sequence  of  images  used  for  case  selection.  We 
further  restricted  the  analysis  to  those  forecasts  where  a  verification  score 
was  available  for  each  forecast  technique.  This  resulted  in  much  smaller 
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sample  sizes  than  those  used  in  the  method  intercomparisons  of  the  previous 
section:  58  cases  (at  increment  1)  to  31  cases  (at  increment  5),  from  the 
forecasts  tor  all  four  seasons.  The  results  shown  in  Figure  55  confirm  the 
conclusion  of  the  previous  section  that  the  segmentation  technique  is  the  best 
of  the  three  extrapolation  techniques  (it  is  significantly  better  at  the  5%  level 
than  any  other  scheme  only  for  increment  1).  In  addition,  it  is  significantly 
better  than  persistence  (at  the  5%  level  of  statistical  significance)  at  all 
forecast  increments.  Movable  persistence  based  on  the  center  of  mass  is  best 
overall,  and  significantly  better  than  persistence  at  all  forecast  increments. 
The  differences  between  the  two  movable  persistence  methods  and  the 
segmentation  technique  are  not  statistically  significant  at  any  increment 
computed  from  this  sample. 

An  examination  of  the  sample  used  for  the  intercomparison  revealed 
that  the  whole  contour  technique  accounted  for  the  majority  of  missing 
forecast  scores.  The  reason  for  these  missing  scores  is  not  entirely  clear 
without  further  investigation.  Since  the  results  so  far  indicate  none  of  the 
three  contour  techniques  are  candidates  for  technology  transition,  this 
problem  was  not  investigated  further;  instead,  it  was  circumvented  by 
deleting  the  whole  contour  technique  from  the  intercomparison.  This 
resulted  in  roughly  doubled  sample  sizes  (140  cases  of  increment  1  to  71  cases 
at  increment  5),  and  correspondingly  tighter  confidence  intervals.  As  is 
shown  in  Fig.  56,  the  movable  persistence  based  on  the  center  of  mass  has 
significantly  better  CSI  scores  than  the  segmentation  technique  at  all 
forecast  increments  (the  difference  between  the  two  movable  persistence 
scores  is  still  statistically  insignificant).  Interestingly,  its  POD  scores  are 
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lower  than  those  of  segmentation,  but  this  is  more  than  compensated  by 
lower  (better)  FAR  scores. 

Finally,  we  compared  contour  forecasts  produced  by  the  cross¬ 
correlation  technique  with  those  of  movable  persistence  (again  based  on  the 
center  of  mass)  for  40  cases  (10  from  each  season).  To  simplify  the 
computations,  we  used  a  slightly  different  verification  procedure  in  this 
comparison  than  in  the  other  contour  forecasts.  In  cases  where  more  than 
one  contour  of  the  selected  contour  level  exists  in  the  scene,  we  treat  all  pixels 
enclosed  by  those  contours  as  cloudy,  and  all  others  as  clear;  previously,  only 
pixels  within  the  one  contour  selected  for  forecasting  were  considered  cloudy. 
For  the  movable  persistence,  the  centroid  position  is  still  defined  in  terms  of 
the  contour  selected  for  forecasting.  The  results  from  this  comparison  (Figure 
57)  show  that  the  cross-correlation  forecasts  have  equal  or  greater  skill  than 
movable  persistence  at  all  forecast  increments.  Out  to  one  hour,  they  are 
significantly  better. 

4.  Summary  of  results  and  recommendations 

The  exhaustive  testing  of  the  contour  extrapolation  techniques  resulted 
in  several  important  findings: 

1)  Identifying,  tracking,  and  forecasting  features  of  interest  on  satellite 
imagery  through  contours  of  brightness  temperature  is  successful  in  only  a 
minority  of  cases.  On  average,  contours  selected  for  forecasting  either  split, 
merged,  or  dissipated  after  three  to  four  hours.  Contours  remained  intact  to 
the  end  of  the  six  hour  sequence  of  satellite  images  used  in  our  study  in  only 
8%  to  33%  of  the  cases,  depending  on  the  season. 
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2)  There  is  no  measurable  skill  in  predicting  anything  but  the  position 
of  contours.  Attempts  at  predicting  smaller  scale  features  (i.e.,  shape)  did  not 
improve  verification  scores,  and  resulted  in  unrealistic  forecast  contour 
shapes.  In  the  case  of  the  segmentation  and  Kawas  techniques,  adjacent 
segment  line  lengths  are  predicted  independently  from  one  another,  resulting 
in  forecast  contours  with  long,  pointed  protrusions.1  In  the  case  of  the  whole 
contour  technique,  the  most  skillful  forecasts  (in  terms  of  verification  scores) 
were  those  with  ellipsoidal  contours. 

3)  Consistent  with  point  2)  above,  the  most  skillful  contour 
extrapolation  technique  is  movable  persistence,  in  which  the  contour  shape  is 
persisted,  and  the  centroid  position  (based  on  the  center  of  mass  of  the  area 
enclosed  by  the  contour)  is  extrapolated  based  on  its  past  movement.  Its 
verification  scores,  along  with  those  of  the  segmentation  technique,  are 
significantly  better  than  persistence  for  forecasts  from  1/2  to  2  1/2  hours  in 
advance. 

4)  The  most  generally  applicable  technique  for  producing  short-term 
forecasts  from  satellite  imagery  is  the  cross-correlation  technique  developed 
in  this  project.  Pixel-based  verification  of  the  cross-correlation  forecasts  show 
an  improvement  against  persistence  out  to  2  1/2  hours  (the  longest  time 
period  studied  for  this  technique).  Contour  forecasts  derived  from  the  cross¬ 
correlation  were  significantly  more  skillful  thar.  those  of  the  best  contour 


1Because  these  protrusions  cover  only  a  small  area,  these  unrealistic  contour  shapes 
are  not  always  accompanied  by  lower  objective  skill  scores. 


33 


extrapolation  technique  (movable  persistence)  out  to  1  hour,  and  no  worse  out 
to  2  1/2  hours. 

The  above  findings  could  possibly  be  altered  if  algorithmic  improve¬ 
ments  to  the  contour  extrapolation  were  explored  (e.g.,  considering  segment 
line  lengths  together,  not  independent  from  one  another).  However,  because 
of  the  overall  failure  rate  of  contour  techniques,  it  is  recommended  that  the 
cross-correlation  technique  developed  in  this  project  be  considered  for 
technology  transition  to  AWDS.  It  requires  only  limited  data  (two  frames  of 
IR  satellite  imagery  spaced  1/2  horn  apart)  and  modest  computer  resources. 

It  will  produce  a  forecast  satellite  image  under  all  conditions,  and  requires 
minimal  human  intervention.  Its  skill  has  been  demonstrated  by  comparison 
to  simple  persistence  in  terms  of  pixel  grayshade  values,  and  against  the  best 
contour  extrapolation  scheme  in  terms  of  contour  forecasts.  It  was  found  to 
yield  the  greatest  improvement  over  persistence  in  cases  dominated  by 
advection;  in  convective  situations,  its  inability  to  handle  development  and 
dissipation  led  to  smaller  improvements  over  persistence.  For  situations  in 
which  contour-based  forecasts  are  useful  (well-defined,  long-lived,  isolated 
features),  and  when  the  additional  forecaster  time  required  for  selection  of 
contour  levels  is  acceptable,  one  may  in  addition  consider  including  the 
movable  persistence  contour  extrapolation  scheme  as  a  candidate  for 
transition.  Although  its  skill  was  found  to  be  inferior  to  that  of  the  cross¬ 
correlation  scheme,  its  computational  simplicity  makes  it  an  attractive 
option. 
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Figure  2:  Domain  used  in  testing  the  cloud  forecast  techniques.  The  domain 
is  a  Lambert-conformal  projection  true  at  30  and  60  N. 


Figure  3:  A  sample  frame  of  raw  GOES  IR  imagery  remapped  to  the  Lambert 
conformal  projection.  Map  overlay  is  not  included  here,  but  the 
domain  corresponds  to  the  one  shown  in  Figure  2.  The  imagery  was 
from  1500  UTC  20  Nov  1991. 
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Figure  5:  The  same  subframe  of  GOES  imagery  as  shown  in  Figure  4,  but 
after  the  selection  of  the  one  feature  for  tracking  and  deletion  of  all 
other  features. 


Figure  4:  A  subframe  of  GOES  imagery  taken  from  3,  with  contours  overlaid 
for  pixel  brightness  181. 
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Stan  Point 


Figure  7:  Illustration  of  the  concept  of  defining  a  contour  using  a  centroid  and 
a  set  of  line  lengths  radiating  from  the  centroid. 
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Figure  8:  Illustration  of  correlation  analysis  and  derivation  of  the 

displacement  vector.  The  displacement  vector  is  directed  from  the 
center  of  the  correlation  box  in  Image  1  to  the  center  of  the  correlation 
box  in  Image  2  most  highly  correlated  with  Image  1. 
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Figure  9:  Example  of  a  correlation  field  used  in  the  derivation  of  a 

displacement  vector.  The  point  (12,2)  is  the  center  of  the  original 
correlation  box.  The  displacement  vector  originates  at  this  point  and 
ends  at  the  point  with  the  highest  correlation. 
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Figure  10:  (a)  Raw  displacement  vectors,  (b)  displacement  vectors  after  a 
quality  control,  and  (c)  objectively  analyzed  displacement  vectors  for 
the  same  satellite  scene  as  shown  in  Figure  5. 
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Figure  11:  Illustration  of  two  contours  bounding  a  forecast  area  (A+C),  and  a 
verification  area  (A  +  B),  and  the  breakdown  into  regions  of  hits  (A), 
misses  (B),  and  false  alarms  (C). 
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Figure  12:  Verification  scores  for  the  CONTROL  sample  for  fall,  for  the  whole 
contour  method  with  method  number  =  2.  Shown  are  CSI  (a),  POD 
(b),  FAR  (c),  BIAS  (d),  and  RMSE  (e)  as  a  function  of  frequency 
(number  of  modes).  The  different  lines  correspond  to  different  forecast 
increments:  increment  1  (solid),  increment  2  (dotted),  increment  3 
(short  dashes),  increment  4  (dash-dot),  increment  5  (dash-dot-dot-dot), 
and  increment  6  (long  dashes). 
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Figure  12:  Verification  scores  for  the  CONTROL  sample  for  fall,  for  the  whole 
contour  method  with  method  number  =  2.  Shown  are  CSI  (a),  POD 
(b),  FAR  (c),  BIAS  (d),  and  RMSE  (e)  as  a  function  of  frequency 
(number  of  modes).  The  different  lines  correspond  to  different  forecast 
increments:  increment  1  (solid),  increment  2  (dotted),  increment  3 
(short  dashes),  increment  4  (dash-dot),  increment  5  (dash-dot-dot-dot), 
and  increment  6  Gong  dashes). 
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Figure  13:  Same  as  Figure  12 


except  for  method  number  =  3 
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Figure  13:  Same  as  Figure  12,  except  for  method  number  =  3. 
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Figure  14:  Same  as  Figure  12,  except  for  method  number  =  4. 
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Figure  15:  Verification  scores  for  the  SLOW  feature  sample  for  fall,  for  the 

method  with  method  number  =  2.  Shown  are  CSI  (a) 
and  RMSE  (b)  as  a  function  of  frequency  (number  of  modes).  The 
different  lines  correspond  to  different  forecast  increments,  as  in 

rigure  12. 
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Figure  16:  Same  as  Figure  15,  except  for  the  FAST  feature  subsample 
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Figure  17:  Same  as  Figure  15,  except  for  the  LARGE  feature  subsample. 


_  Inc  1 

.  Inc  2 

- Inc  3 

_  Inc  4 

-  Inc  5 

- Inc  6 


Number  of  Modes 


b 

Whole  Contour  PRE2  FALL  Method  No.=  2  RMSE 


Figure  19:  Same  as  Figure  15 


except  for  the  PKE2  subsampl, 


Figure  21:  Same  as  Figure  15,  except  for  the  PRE6  subsample, 
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Figure  22:  Same  as  Figure  12,  except  for  winter. 
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Figure  23:  Same  as  Figure  12,  except  for  spring. 
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Figure  24.  Same  as  Figure  12,  except  for  summer. 
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Figure  24:  Same  as  Figure  12,  except  for  summer. 
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Figure  25:  Verification  scores  for  the  CONTROL  sample  for  fall,  for  the 

S-h0Wn  are  CSI  (a)’  P0D  «>).  FAR  (c),  BIAS  (d), 
and  RMSE  (e)  as  a  function  of  the  segmentation  angle  (in  °)  The 

different  lines  correspond  to  different  forecast  increments,  as  in 
rigure  12. 
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Figure  25:  Verification  scores  for  the  CONTROL  sample  for  fall,  for  the 

segmentation  method.  Shown  are  CSI  (a),  POD  (b),  FAR  (c),  BIAS  (d), 
and  RMSE  (e)  as  a  function  of  the  segmentation  angle  (in  °).  The 
different  lines  correspond  to  different  forecast  increments,  as  in 
Figure  12. 
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Figure  26:  Verification  scores  for  the  SLOW  feature  sample  for  fall,  for  the 

segmentation  method.  Shown  are  CSI  (a)  and  RMSE  (b)  as  a  function 
of  the  segmentation  angle  (in  °).  The  different  lines  correspond  to 
different  forecast  increments,  as  in  Figure  12. 
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Figure  30:  Same  as  Figure  26,  except  for  the  PRE2  subsample. 
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Figure  32:  Same  as  Figure 
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Figure  33:  Same  as  Figure  25,  except  for  winter. 
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Figure  34:  Same  as  Figure  25,  except  for  spring. 
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FALL  KAVVAS  CONTROL  CSI 
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b 

FALL  KAVVAS  CONTROL  POD 
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Figure  36:  Verification  scores  for  the  CONTROL  sample  for  fall,  for  the 

Kawas  method.  Shown  are  CSI  (a),  POD  (b),  FAR  (c),  BIAS  (d),  and 
RMSE  (e)  as  a  function  of  the  smoothing  coefficient.  The  different 
lines  correspond  to  different  forecast  increments,  as  in  Figure  12. 
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Figure  36:  Verification  scores  for  the  CONTROL  sample  for  fall,  for  the 

Kawas  method.  Shown  are  CSI  (a),  POD  (b),  FAR  (c),  BIAS  (d),  and 
RMSE  (e)  as  a  function  of  the  smoothing  coefficient.  The  different 
lines  correspond  to  different  forecast  increments,  as  in  Figure  12. 
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Figure  37:  Verification  scores  for  the  SLOW  feature  sample  for  fall,  for  the 

Kawas  method.  Shown  are  CSI  (a)  and  RMSE  (b)  as  a  function  of  the 
smoothing  coefficient.  The  different  lines  correspond  to  different 
forecast  increments,  as  in  Figure  12. 
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FALL  KAVVAS  FAST  CSI 


Figure  38:  Same  as  Figure  37,  except  for  the  FAST  feature  subsample. 
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Figure  39:  Same  as  Figure  37,  except  for  the  LARGE  feature  subsample 
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Figure  41:  Same  as  Figure  37,  except  for  the  PRE2  subsample. 
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Figure  42:  Same  as  Figure  37,  except  for  the  PRE4  subsample 


Figure  43:  Same  as  Figure  37,  except  for  the  PRE6  subsample. 
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Figure  44:  Same  as  Figure  36,  except  for  winter. 
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Figure  45:  Same  as  Figure  36,  except  for  spring. 
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Figure  46:  Same  as  Figure  36,  except  for  summer. 
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Figure  46:  Same  as  Figure  36,  except  for  summer. 
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Figure  47:  Verification  scores  for  the  CONTROL  sample  for  fall  as  a  function 

Shown  are  the  sample  size  (a),  CSI  (b)  POD  (c) 
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Figure  47:  Verification  scores  for  the  CONTROL  sample  for  fall  as  a  function 
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and  RMSE  (c) .  The  different  lines  correspond  to  different  forecast 
techniques  ,  as  in  Figure  47 . 


Figure  50:  Same  as  Figure  48,  except  for  the  LARGE  feature  subsample 
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Figure  52:  Same  as  Figure  47,  except  for  winter 
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Figure  with  6%  confidence  intervals  from  the 

AMU  V  A  lor  the  CONTROL  sample  for  the  whole  year  as  a  function  of 
forecast  increment.  Shown  are  the  CSI  (a),  POD  (b),  and  FAR  (c)  The 
different  lines  .and  symbols  used  for  the  confidence interval 
correspond  to  different  forecast  techniques:  segmentation  (solid 

X)’ Kavvas  (dashed.  diamondXand 
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Figure  56:  As  Fig.  55,  but  with  the  whole  contour  technique  eliminated  from 
the  comparison. 
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Mean  FAR  with  5%  confidence  intervals 


Fig.  57:  Mean  verification  scores  with  5%  confidence  intervals  for  the  contour 
forecasts  from  the  cross-correlation  technique  (solid,  asterisk)  and 
movable  persistence  based  on  the  center  of  mass  (dashed,  diamond). 
From  the  ANOVA  for  the  40  case  sample  for  the  whole  year.  Shown 
are  the  CSI  (a),  POD  (b),  and  FAR  (c)  as  a  function  of  forecast 
increment. 
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Appendix  A:  A  SHORT  -  TERM  CLOUD  FORECAST  SCHEME 
USING  CROSS  CORRELATIONS 

Thomas  M.  Hamill  and  Thomas  Nehrkom 


Atmospheric  and  Environmental  Research,  Inc. 
840  Memorial  Dr 
Cambridge,  Massachusetts  02139 


1.  INTRODUCTION 

This  paper  describes  a  trajectory-based  cloud 
forecast  technique  based  on  lag  cross  correlations 
(Panofsky  and  Brier  1968).  This  technique  generates 
a  set  of  synthetic  IK  images  nearly  indistinguishable 
from  the  real  images.  It  is  suitable  for  short-range 
forecasts  and  can  be  run  in  a  few  minutes  on  current 
generation  workstations.  This  scheme  requires  two 
IR  satellite  images  of  the  same  scene  not  more  than 
one  half  hour  apart;  no  other  data  are  needed. 
Advective  velocities  are  generated  through  the  cross¬ 
correlation  of  multiple  subsets  of  the  two  frames. 
This  is  based  on  the  same  technique  used  for  satellite- 
derived  winds  (Merrill  1989),  and  has  also  been  used 
for  cloud  and  rainfall  forecasts  by  Bellon  et  al.  (1992) 
and  applied  to  radar  for  boundary  layer  winds  by 
Tuttle  and  Foote  (1990).  This  particular  application 
was  designed  with  Air  Force's  Automated  Weather 
Distribution  System  workstation  in  mind,  but  should 
be  similarly  useful  within  the  civilian  community. 

Since  the  technique  described  here  is  trajectory- 
based,  it  has  some  typical  problems  associated  with 
such  techniques.  For  example,  it  cannot  effectively 
develop  or  dissipate  clouds,  limiting  the  length  of 
time  and  the  weather  regime  for  which  they  are 
useful.  Another  problem  is  determining  correct 
advecting  velocities;  if  an  arbitrary  level  (say,  500 
mb)  wind  field  is  used,  low  clouds  may  move  unduly 
fast  and  high  clouds  too  slowly.  As  will  be  shown,  a 
correlations-based  approach  partly  avoids  this  latter 
problem. 

Longer-range  forecast  schemes  (e.g.,  Slingo, 
1987;  Sundqvist  et  al.  1989;  Mitchell  and  Hahn 
1989)  typically  require  the  execution  of  a  dynamic 
forecast  model.  They  have  their  own  set  of 
drawbacks,  such  as  poor  spatial  resolution,  sensitivity 
to  model  spinup,  and  dependence  on  hard-to-forecast 
variables  such  as  relative  humidity.  There  is  a 
crossover  point  between  6  and  18  hours  where  these 
longer-range  forecast  schemes  become  more  skillful 
than  trajectory  schemes.  Our  technique  was  designed 
to  meet  the  very  short-range  forecast  (0-2  1/2  h) 
needs  of  the  weather  consumer,  as  such,  it  is  clearly 
more  skillful  than  any  of  the  longer-range  schemes 
during  this  period.  Our  testing  was  limited  to  this 


forecast  envelope;  thus,  we  have  not  yet  developed 
any  statistical  basis  for  determining  the  crossover 
time  when  this  scheme's  skill  is  exceeded  by  one  of 
the  longer-range  schemes. 

2.  ALGORITHM  DESCRIPTION 

2.1  Review  of  Cross-Correlations  Technique 

At  the  heart  of  our  prototype  cloud  forecast 
technique  is  the  derivation  of  displacements  vectors 
(i.e.,  "winds")  through  a  cross-correlations  analysis. 
The  technique  is  as  follows;  a  subset  of  pixels  from 
the  First  image  in  a  satellite  loop  is  chosen.  For 
purpose  of  illustration  in  Figure  1,  this  subset  is  8*8; 
in  our  actual  prototype  scheme,  the  subset  chosen  is 
15*15.  Next,  based  on  the  maximum  possible  wind 
displacement  in  half  an  hour,  a  search  radius  is 
chosen,  and  identically  sized  subsets  of  pixels  from 
the  second  image  with  centers  inside  the  search 
radius  are  correlated  against  the  subset  from  the  first 
image.  If  the  scene  has  a  uniform  cloud  layer,  then 
the  advective  velocity  for  this  layer  will  be  defined  by 
the  vector  from  the  center  of  the  subset  in  Image  1  to 
the  center  of  the  subset  with  the  maximum  correlation 
in  Image  2.  This  is  illustrated  in  Figure  1,  and  a 
sample  plot  of  correlation  coefFicients  and  a  derived 
displacement  vector  is  shown  in  Figure  2.  In  our 
technique,  a  cross  correlations  analysis  is  repeated  at 
a  gridded  subset  of  points  throughout  the  domain,  and 
wind  displacement  vectors  for  all  points  are  then 
derived  through  an  objective  analysis.  {The  domain 
is  a  lambert-conformal  map  projection  covering  the 
eastern  two-thirds  of  the  U.S.  and  southern  Canada. 
Raw  GOES  IR  imagery  was  remapped  into  a 
256*256  set  of  pixels  in  this  projection,  with  an 
approximate  resolution  of  11  km.] 

The  correlations  technique  has  drawbacks;  there 
may  be  multiple  cloud  layers,  in  which  case  a  single 
displacement  vector  may  not  be  applicable  for  that 
region  (Leese  and  Novak  1971).  Second,  if  the  frame 
is  clear,  then  the  displacement  vector  will  be  null, 
even  with  strong  winds;  use  of  the  null  wind  in  a 
trajectory  analysis  will  prevent  clouds  from  moving 
into  the  clear  region.  Further,  if  small  subsets  are 
used  for  the  correlation  analysis,  the  displacement 
vector  may  be  inaccurate,  reflecting  a  random  high 
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correlation  and  not  a  true  advectivc  wind  velocity. 

The  focus  here  is  on  producing  a  cloud  forecast 
and  not  on  the  derivation  of  accurate  wind 
observations  for  use  in  forecast  i  .wdel  initialization; 
thus,  these  limitations  can  either  be  corrected  or 
overcome.  For  the  first  problem,  a  displacement 
vector  for  a  region  with  multiple  cloud  layers  may 
represent  die  most  appropriate  single  displacement 
vector  that  can  be  used  in  a  one-layer  trajectory 
scheme.  Further,  let  us  assume  advcction  (rather  than 
development  or  dissipation)  is  the  dominant  process, 
and  there  is  only  one  cloud  layer  in  a  given  area.  If 
so,  then  this  technique  will  derive  more  reasonable 

SEARCH  RADtUX  BOX  WITH  MAX 


Figure  1.  Illustration  of  correlation  analysis  «nd  derivation 
of  the  displacement  vector.  The  displacement  vector  is 
directed  from  the  center  of  the  correlation  box  in  Image  1  to 
the  center  of  the  correlation  box  in  Image  2  most  highly 
correlated  with  Image  1.  . 


0  4  8  12  16  20  24 

x-pixel  location  (I2=center) 


Figure  2.  Example  of  a  correlation  field  used  in  (he 
derivation  of  a  displacement  vector.  The  point  (122)  is  the 
center  of  the  original  correlation  box.  The  displacement 
vector  originates  at  this  point  and  ends  at  the  point  with  the 
highest  correlation. 


vectors  than  could  be  achieved  using  just  single  level 
winds  such  as  a  500  mb  wind  field.  Wind 
magnitudes  generally  increase  with  height,  so  by 
using  the  cross -correlations  technique,  areas  of  low 
cloud  will  have  smaller  magnitude  displacement 
vectors  than  neighboring  areas  with  high  cloud.  In 
essence,  the  trajectory  wind  field  over  the  domain  can 
be  warped  to  represent  correct  advectivc  velocities 
for  a  given  region's  cloud  height. 

Should  part  of  the  frame  be  cloud-free,  those 
areas  will  have  null  displacement  vectors.  Since  the 
cloud  edge  moving  into  an  area  is  often  one  of  the 
most  crucial  features  to  forecast,  our  technique  was 
tailored  to  avoid  the  anomalous  slowing  of  cloud 
edges.  Hence,  a  restriction  was  imposed  during  the 
objective  analysis  step,  where  vectors  are  set  for  all 
pixels  in  the  domain.  The  imposed  restriction 
removes  null  displacement  vectors  from  use  as 
observations.  As  a  result,  wind  velocities  at  cloud 
edges  are  not  decreased  from  the  assimilation  of  the 
null  displacement  vectors,  and  cloud  edges  are 
advected  at  more  appropriate  speeds.  However, 
without  correction,  clear  areas  in  front  of  the 
advancing  clouds  are.also  advected.  This  presents  no 
problem  when  the  clear-scene  background  is 
homogeneous,  but  if  terrain  features  or  lakes  are 
present,  they  may  appear  to  move.  This  problem 
was  alleviated  by  identifying  dear  and  doudy  pixels 
(using  user-defined  threshold  values)  and  using 
persistence  rather  than  advcction  for  dear  areas. 

The  last  problem,  inaccurate  vectors  and 
resultant  displacement  field  noisiness,  can  be 
corrected  in  two  ways:  first,  the  size  of  the  subset  of 
pixels  for  correlation  can  be  increased.  However,  this 
is  done  at  the  ride  Of  increased  computations.  Our 
choice  of  a  15*15  subset  was  a  compromise  between 
excessive  CPU  time  and  excessive  noise.  The  second 
correction  used  here  is  a  local  consistency  check  to 
compare  each  derived  vector  against  its  neighbors.  In 
this  scheme,  displacement  vectors  that  deviate  more 
than  the  width  of  four  pixels  from  the  local  mean 
vector  are  replaced  by  the  local  mean  (other 
deviations  were  tried,  but  four  pixels  gave  the  best 
result  in  our  tests).  Figure  3  shows  a  frame  of  the 
remapped  GOES  IR  satellite  imagery  used  in  the 
correlation  (1230  UTC  20  November  1991).  Figure  4 
shows  a  field  of  correlation  displacement  vectors 
derived  from  data  between  1200  and  1230  UTC  20 
November  1991,  illustrating  typical  noise.  Figure  5 
shows  the  same  field  of  displacement  vectors  after 
application  of  the.  consistency  check.  This  set  of 
vectors  is  itself  quite  consistent;  for  example, 
compare  the  results  in  Figure  5  with  those  in  Figure 
6.  which  shows  the  derived  vectors  for  the  frames  1/2 
hour  later,  from  1230  and  1300  UTC  imageryr 
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The  next  step  after  derivation  and  QC  of  the 
displacement  vectors  is  to  take  the  field  of  vectors 
and  produce  a  gridded  displacement  field  for  every 
point  in  the  domain.  This  is  done  by  objective 
analysis.  The  objective  analysis  technique  and  radius 
of  influence  was  chosen  to  smooth  the  vector  field 
slightly  and  produce  a  continuous  flow  pattern.  We 
used  two-pass  objective  Cressman-type  objective 
analysis  with  a  first  pass  summation  function 
(Benjamin  and  Seaman  1985)  of  the  form: 

nobs 

X  (W  2  *  DCxn.yn» 

Gl(ij)  =  n~I _ 

nobs 

X  (W) 
n=l 


Here  G1  is  the  first  guess  at  pixel  (ij)»  W  is  the 
standard  isotropic  Cresstnan  weighting  function 
dependent  on  the  distance  between  the  observation 
and  the  analysis  point,  and  D  is  the  displacement 
vector  velocity  for  the  observation  at  location  (xn,yn)- 
U  and  V  components  are  analyzed  separately  in  this 
process.  The  second  pass  summation  function  is  of 
the  form: 

nobs 

X  (W2*(D(xn.yn)-Gl(xn,yn)) 

G20J)=  n=] _ 

nobs 

X  (W) 

n=I 

where  G2  represents  the  final  analysis  value.'  The 
radius  of  influence  was  set  at  90  and  40  pixels  for  the 
first  and  second  pass,  respectively.  As  noted  earlier, 
we  throw  out  the  null  vectors  before  proceeding  with 
the  objective  analysis.  Figures  7  and  8  illustrate  the 
effect  of  this;  Figure  7  shows  an  objective  analysis  of 
the  data  in  Figure  5,  allowing  the  null  vectors  to  be 
used.  Figure  8  shows  the  implemented  version,  with 
the  null  vectors  deleted.  As  can  be  seen  in  the  upper- 
left  comer  of  the  domain,  the  influence  of  the 
windfield  is  spread  into  the  clear  areas,  resulting  in 
stronger  magnitude  winds  along  cloud  edges. 


2-3  Semi-La grangian  Trajectory  Forecast 

In  order  to  produce  forecast  images,  pixel 
intensities  from  the  latest  satellite  image  are  now 
advccted  using  the  full  field  of  displacement  vectors. 


0A  Vectors  w/o  Null  Elimination 


Figure  7.  Objective  analysis  of  data  from  Figure  5.  without 
null  displacement  vectors  eliminated. 

0A  Vectors  with  Null  Elimination 


Figure  8.  Objective  analysis  of  data  from  Figure  5,  with 
null  displacement  vectors  eliminated. 


Before  the  trajectory  technique  is  performed, 
however,  the  U  ar.d  V -component  1/2  hour 
displacement  vectors  are  modified  slightly  through  a 
semi-lag  rangian  scheme  (Stani forth  and  Cote  1991). 
The  use  of  semi-lagrangian  displacement  vector 
rather  than  the  original  displacement  vector  itself 
markedly  improves  the  trajectory  forecasts  in  areas 
where  the  wind  field  at  the  trajectory  endpoint  differs 
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significantly  from  the  wind  field  at  the  trajectory 
origin.  In  essence,  this  scheme  iterates  to  find  the 
correct  compromise  trajectory. 

With  final  advective  velocities  calclated,  a 
simple  backwards-in-time  trajectory  is  used  with  the 
initial  satellite  data  to  make  forecast  images.  In  this 
scheme,  the  trajectory  endpoints  for  time  interval  t(i) 
are  known;  they  are  the  regular  set  of  gridpoints  for 
the  domain.  Using  the  semi-lagrangian  vectors,  the 
scheme  finds  the  trajectory  origin  at  time  t(i-l);  if  the 
origin  is  co'  cated  with  a  pixel,  that  pixel  value  is 
used  as  the  forecast  value;  however,  if  the  origin 
does  not  lie  directly  on  a  pixel,  then  a  bilinear 
interpolation  of  the  surrounding  four  pixel  values  is 
performed  to  determine  the  forecast  pixel  value.  If 
the  trajectory  origin  is  outside  the  domain,  the 
trajectory  is  assumed  to  originate  at  the  nearest 
border  point.  The  forecast  scheme  steps  forward  in 
1/2  hour  intervals,  with  each  successive  forecast 
frame  used  in  the  initialization  of  the  next  step. 


3.  RESULTS 

Our  prototype  cloud  forecast  scheme  has  now 
been  tested  with  satellite  images  over  the  central  U.S. 
for  each  season.  It  shows  skill  over  persistence  for 
all  time  periods  tested  (1/2  to  2  1/2  hours  in  1/2  hour 
increments).  The  scheme  was  also  tested  side-by- 
side  with  an  advection  Lsheme  using  the  500  mb  wind 
field,  and  comparison  tests  were  made  against  several 
contour  extrapolation  methods. 


3.1  Comparison  with  Persistence  and  5QQ  mb 

Advection  Forecasts 

For  all  tests  and  all  times,  the  correlations-based 
advection  scheme  showed  less  error  than  persistence. 
An  ensemble  of  10  forecasts  from  each  of  the  four 
seasons  were  made,  and  RMS  errors  (in  pixel 
grayshade  values)  were  calculated  using  the  observed 
satellite  imagery  for  verification.  The  fall,  winter, 
spring,  summer  error  statistics  for  correlations 
forecasts  and  persistence  are  summarized  in  Table  1. 
As  shown,  both  persistence  and  correlations  exhibit 
higher  errors  in  summer  than  in  winter.  From  visual 
inspection,  this  was  clearly  due  to  the  dominance  of 
convective  development  and  dissipation  over 
advection  during  the  summer  months.  Conversely, 
during  the  winter,  advective  processes  seemed  to 
dominate,  and  both  persistence  and  the  correlations 
technique  show  smaller  errors,  though  the 
percentage  improvement  of  correlations-based 
forecasts  over  persistence  is  greater  in  winter  than  in 
summer. 


Fcst. 

Season 

Correlation 

500  mb 

Persistence 

Hour 

RMS 

RMS 

RMS 

1/2 

Fall 

9.13 

— 

1132 

Winter 

8.57 

936 

11.30 

Spring 

8.98 

- 

11.95 

Summer 

11.28 

11.63 

13.96 

1 

Fall 

12.87 

16.21 

Winter 

11.67 

12.82 

14.98 

Spring 

12.55 

- 

16.04 

Summer 

1632 

16.99 

19.77 

11/2 

Fall 

1535 

— 

18.59 

Winter 

14.19 

15.22 

1738 

Spring 

14.99 

- 

1835 

Summer 

19.94 

20.61 

23.46 

2 

Fail 

17.27 

— 

20.48 

Winter 

16.12 

17.09 

19.23 

Spring 

17.02 

- 

2031 

Summer 

22.77 

23.44 

26.36 

21/2 

Fall 

18.95 

— 

22.16 

Winter 

17.79 

18.64 

20.83 

Spring 

18.64 

— 

22.19 

Summer 

25.25 

25.82 

28.78 

Table  1.  Comparison  of  RMS  Errors  for  for  10-case 
ensembles  of  correlation.  500  mb,  and  persistence  forecasts 
over  each  of  the  four  seasons,  and  for  each  forecast  interval. 
Units  are  grayshade  values,  which  range  from  0-256  in  the 
original  imagery. 


A  reasonable  alternative  candidate  to  a 
correlations-based  forecast  scheme  is  one  using  500 
mb  winds  for  displacement  vectors.  To  test  this,  we 
also  compared  a  set  of  500-mb- based  forecasts  to  the 
observed  IR  pixel  values  fra-  the  winter  and  summer 
cases.  The  results  for  these  are  also  presented  in 
Table  1.  Generally,  the  correlations  forecasts  are 
clearly  superior  in  the  winter,  but  much  closer  to  the 
RMS  of  the  500  mb  forecasts  in  the  summer.  The 
generally  higher  RMS  in  the  summer  and  the  lesser 
RMS  difference  between  correlations  and  500  mb 
forecasts  for  summer  cases  is  most  likely  due  to  the 
dominance  of  convective  development  and 
dissipation.  During  the  summer,  use  of  the  500  mb 
wind  field  as  a  steering  current  for  thunderstorms 
generally  makes  as  much  sense  as  correlations- 
derived  winds,  which  cannot  be  expected  to  produce 
reliable  wind  vectors  in  regions  of  development  and 
dissipation.  Conversely,  during  the  winter  months, 
the  correlation  scheme's  ability  to  determine  an 
accurate  advecting  velocity  regardless  of  the  height  of 
the  cloud  is  likely  responsible  for  its  higher  skill. 

A  few  forecasts  were  made  using  700  mb  winds 
instead  of  500  mb;  the  verification  scores  for  this 
small  sample  were  generally  worse  using  700  mb 
data  than  500  mb. 
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Other  techniques  exist  for  the  extrapolation  of 
prominent  cloud  features,  notably  contour 
extrapolation  methods  (Heideman  et  al.  1990; 
Ruggiero  et  al.  1991).  For  these  methods,  significant 
weather  features  such  as  thunderstorms  are  contoured 
at  a  user-selected  brightness  level,  and  the  past 
behavior  of  the  contour  shape  and  movement  is  used 
to  forecast  future  shape  and  movement.  The  cross¬ 
correlations  based  forecasts  were  tailored  to  also 
produce  contour  forecasts  and  are  currently  being 
compared  side-by-side  against  the  existing  contour 
forecast  schemes.  The  critical  success  index,  false 
alarm  ratio,  and  probability  of  detection  (Stanslci  et 
al.  1989)  will  be  the  benchmark  performance  indices 
in  this  comparison.  Early  results  indicate  competitive 
performance  compared  with  contour-based  schemes. 
Quantitative  results  will  be  available  at  the 
conference. 


4.  CONCLUSIONS 

A  cross-correlations  based  cloud  forecast 
technique  has  been  developed  which  is  skillful  at 
extrapolating  cloud  features.  The  output  is  a  set  of 
synthetic  satellite  images  which  may  be  looped  in 
combination  with  analyzed  images.  This  scheme  was 
designed  for  a  U.S.  Air  Force  weather  workstation, 
but  it  may  also  be  useful  to  forecasters  and  broadcast 
meteorologists  who  need  short-range  cloud  forecasts. 
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